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Gas driven massive black hole binaries: signatures in the 
nHz gravitational wave background 
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ABSTRACT 

Pulsar timing arrays (PTAs) measure nHz frequency gravitational waves (GWs) gen- 
erated by orbiting massive black hole binaries (MBHBs) with periods between 0.1-10 
yr. Previous studies on the nHz GW background assumed that the inspiral is purely 
driven by GWs. However, torques generated by a gaseous disk can shrink the binary 
much more efficiently than GW emission, reducing the number of binaries at these sep- 
arations. We use simple disk models for the circumbinary gas and for the binary-disk 
interaction to follow the orbital decay of MBHBs through physically distinct regions of 
the disk, until GWs take over their evolution. We extract MBHB cosmological merger 
rates from the Millennium simulation, generate Monte Carlo realizations of a popu- 
lation of gas driven binaries, and calculate the corresponding GW amplitudes of the 
most luminous individual binaries and the stochastic GW background. For steady- 
state ev-disks with a > 0.1 we find that the nHz GW background can be significantly 
modified. The number of resolvable binaries is however not changed by the presence 
of gas; we predict 1-10 individually resolvable sources to stand above the noise for 
a 1-50 ns timing precision. Gas driven migration reduces predominantly the number 
of small total mass or unequal mass ratio binaries, which leads to the attenuation 
of the mean stochastic GW-background, but increases the detection significance of 
individually resolvable binaries. The results are sensitive to the model of binary-disk 
interactio n. The GW background is not attenua ted significantly for time-dependent 
models of llvanov. Papaloizou. fc Polnarevl (jl999t ). 
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1 INTRODUCTION 



Inspiraling massive black hole binaries (MBHBs) with 
masses in the range ~ 10 4 — 10 10 Mq are expected to be the 
dominant source of gravitational waves (GWs) at ~ nHz 
- mHz frequencies ilHaehnelt I [l99i; Ijaffe fc Backer 1 120031 ; 
IWvithe fc Loeb 1 120031 ; ISesana et all |2004 l2005h . The fre- 
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The com- 
the Euro- 



the Parkes radio-telescope 
plete Parkes PTA (PPTA _ 
pean Pulsa r Timing Array (EP TA, Janssen et alj200Sl ). and 
NanoGrav jjenet et al. 1120091 ) are expected to improve con- 
siderably on the capabilities of these surveys, eventually 
oining their efforts in the international PTA project (IPTA, 



quency band ~ 10 Hz - 1 Hz will be probed by the Lase r 
Interferometer Space Antenna (LISA, feender et al.l Il998l ). 
a space-borne GW laser interferometer developed by ESA 
and NASA. The observational window 10" 9 Hz - 10" 6 Hz, 
corresponding roughly to orbital periods 0.03 - 30 yr, is al- 
ready accessible using Pulsar Timing Arrays (PTAs; e.g. 
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iHobbs et aL 2010[); an d the planned Square Kilometer Ar 



ray (SKA. lLazio 1 120091 s ) will produce a major leap in sensi- 
tivity. 

Radio pulses generated by rotating neutron stars travel 
through the Galactic interstellar medium and are detected 
by radio telescopes on Earth. The arrival times of pulses 
are fitted for a model including all the known and measured 
systematic effects a ffecting the signal generation, propaga - 
tion and detection (|Edwards. Hobbs. fc Manchester 1120061 s ]. 
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Timing residuals between the observed pulses and the best 
fit model, carry information on additional unmodelled ef- 
fects, including the presence of GWs. Indeed, GWs mod- 
ify the propagation of radio signals from the pulsar to the 
Earth (ISazhin I Il978l: [Petweiler Ill979l; iBertotti et al.lll983l ; 



iHellines fc Downs Ill983l : ljenet et al. l2005h . and PTAs mea- 
sure the direction dependent systematic variations in the 
arrival times of signals from a sample of nearly stationary 
pulsars in the Galaxy distributed over the sky. 

PTAs provide a direct observational window onto the 
MBH binary population, and can contribute to address a 
number of open astrophysical questions, such as the shape 
of the bright end of the MBH mass function, the nature of 
the MBH-bulge relation at high masses, and the dynamical 
evolution at sub-parsec scales of the most massive binaries 
in the Universe (par ticularly relevant to the so-cal led "fi- 
nal parsec problem", iMilosavlievic fc Merritt"ll2003l ) PTAs 
can detect gravitational radiation of two forms: (i) the 
stochastic GW background produced by the incoherent su- 
perposition of radiation from the whole cosmic population 
of MBHBs and (ii) individual sources that are sufficiently 
bright in GWs to outshine the background (typically mas- 
sive, M>1O 9 M0, and "cosmologically nearby", dz,<3Gpc). 
Both classes of signals are of great interest, and PTAs could 
lead to the discovery of systems difficult to detect with other 
techniques (for alternatives using active galactic nuclei, see 
lHaiman. Kocsis. fc Menou 1120091 . and references therein). 

Popular scenarios of ma ssive black hole (MBH) for- 
mation and evolution (e.g. Volonteri. Haardt. fc Madaul 



20031; iKoushiappas fc Zentner I 120061 ; iMalbon et all bool 



Yoo et al. 1120071) predict frequent MBH mergers (up to sev 



eral hundreds per year), implying the existence of a large 
number of sub-parsec MBHBs. The prospect for detecting 
GW signals using PTAs depends on the number and cosmo- 
logical distribution of MBHBs with orbital periods of 0.03 - 
30 yr, or separations typically in the range 0.001 — 0.1 pc. The 
three main ingredients for calculating the GW background 
are 

(i) the merger rate of MBHBs as a function of mass and 
redshift, 

(ii) the relative time each binary spends at these separa- 
tions during a merger episode and 

(iii) the amplitude of the GW signal produced by each 
individual stationary system. 



Recently ISesana. Vecchio. fc Colacino I (|2008| . SVC08 

(|2009l . 



hereinafter) and ISesana. Vecchio. fc Volonteri 



SVV09 hereinafter), carried out a detailed study of the 
expected signals (stochastic and individual) , focusing on the 
uncertainties related to (i). They found that the background 
is affected by the galaxy merger rate evolution along the 
cosmic history, the massive black hole mass function, and 
the accretion history of the MBHB during a galaxy merger, 
and they predict a factor of ~ 10 uncertainty for the charac- 
teristic strain amplitude in the range 2 x 10 -15 — 2 x 10 -16 , 
at / = 1/yr, within the expected detection capabilities of 
the complete PPTA and of the SKA. They pointed out that 
the GW signal can be separated into individually resolvable 
sources and a stochastic background, and found the number 
of individually resolvable sources for a Ins timing precision 
level to be between 5 to 15, depending on the considered 
model. 



In this paper, we examine for the first time how pre- 
dictions relevant for PTA observations are modified by the 
presence of ambient gas, affecting the inspiral rate of bina- 
ries during a merger episode, ingredient (ii) above. A gaseous 
envelope is expected to surround the binary because MB- 
HBs are produced in galaxy mergers, which are known to 
trigger inflows of large quantities of gas into the central re- 
gion, as shown by hydrodynamic simulations (|Springel et alj 
l2005h . This gas, accreted onto the MBHs, is responsible 
for luminous AGN activity, and is also expected to cat- 
alyze the coalescence of the new-for med MBH pair (e.g., 
lEscala et all 12004 iDotti et al]|2007l) . as described below. 
The forming MBHB spirals inward initially as a result of 
dynamical friction on dark matter, ambient stars, and gas 
l|Begelman. Blandford. fc Reeslll980t ) . As the binary separa- 
tion shrinks to sub-pc scales, the supply rate of stars crossing 
the orbit decreases, and the interaction with stars becomes 
less and less efficient to shrink the binary. In gas rich merg- 
ers, the dense nuclear gas is expected to cool rapidly and 
settle into a geomet r ically-thin circumb inary accretion disk 
(e.g., iBarnes 1 12001 lEscala et al"1l2004 ). Torques from the 
tidal field of the binary clear a gap in the gas with a radius 
less than twice the separation of the binary, and generates a 
spiral density wave in the gaseous disk, which in turn drains 
angular momentum away from the binary on a relatively 
short timescale within the last parsec, <10 7 vr (lEscala et alj 
20051; lArmitage fc Nataraian 112001 120051, iDotti et al.ll2007h 



Havasaki II2009I . see however Lodato et al 12009? ). Ultimately, 
at even smaller separations, corresponding to an orbital 
timescale of ~ years, the emission of GWs becomes the dom- 
inant mechanism driving the binary to the final coalescence. 
The main point of this paper, is to notice that the most 
sensitive PTA frequency band corresponds to orbital sepa- 
rations near the transition between gas and GW dominated 
evolution . As binaries shrink more quickly inwards in the 
gas-driven phase, the number of binaries emitting at each 
given separation is decreased compared to the purely GW- 
driven case. The subject of this work is to explore how var- 
ious gas-driven models modify the expectations on the GW 
signal potentially observable by PTAs. 

Recently, lHaiman. Kocsis. fc Menou I (|2009l . HKM09 
hereafter), examined the evolution of MBHBs in the 
gas-driven regime for simple mod els of geometrically 
thin circumbinary disks (se e also, ISver fc Clarke! 1 19951 : 
lArmitage fc Nataraian 1120051 ). The interaction between the 
binary and the gaseous disk is analogous to type-II plane- 
tary migration, and evolves through two main phases. First, 
the inspiral is analogous to the disk-dominated type-II mi- 
gration of planetary dynamics, where the binary migrates 
inwards with a radial velocity equal to that of the gas accret- 
ing towards the center. Later, as the mass of the gas within a 
few binary separations becomes less than the reduced mass 
of the binary, the evolution slows down, and it is analogous 
to the planet-dominated (or secondary-dominated) type- 
II migration. In both cases, the radial inspiral rate is still 
much faster than in the purely GW driven case at orbital 
separations beyond a few hundred Schwarzschild radii. For 
standard Shakura-Sunyaev a-disk models, the viscosity is 
assumed to be proportional to the total pressure, and is 
consequently very large in the radiation pressure dominated 
phase at small radii, increasing the migration rate in the 
radiation pressure dominated phase. On the other hand, for 
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/3-disk models, where the viscosity is proportional to the 
gas pressure only, the increase of radiation pressure does 
not impact the viscous timescale, and the migration rate is 
relatively slower in this regime. Finally, we note that the 
binary-gas interaction is al so significantly different for non- 
stead y models of accretion (jlvanov. Papaloizou. fc Polnarevl 
1 19991 . HKM09). In the typical secondary-dominated phase, 
gas flows in more quickly then how the binary separation 
shrinks, and is repelled close to the outer edge of the gap 
by the torques of the binary. This causes gas to accumulate 
near the gap, and delays the merger of the binary relative 
to the steady-state models. 

In this paper, we couple the HKM09 models for the 
migration of MBHBs in the presence of a steady gaseous 
disk, to the population models derived in SVV09, and we 
compute the effects on GWs at nHz frequencies. Here, 
we restrict to nearly circular inspirals for simplicity, us- 
ing the corresponding GW spectrum (ingredient-iii above). 
This assumption might be violated in gas driven inspirals 
I Armitage fc Nataraian 1 120051 ; ICuadra et al.l l2009h . and is 
the subject of a future paper (Sesana & Kocsis 2010, in 
prep). 

The paper is organized as follows. In Section 2 we in- 
troduce the theory of the GW signal from a MBHB popula- 
tion, describing its characterization in terms of its stochastic 
level and of the statistics of individually resolvable sources. 
In Section 3 we describe our MBHB population model, cou- 
pling models of coalescing binaries derived by cosmological 
N-body simulations to a scheme for the dynamical evolution 
of the binaries in massive circumbinary disks. We present in 
detail our results in Section 4, and in Section 5 we briefly 
summarize our main findings. Throughout the paper we use 
geometric units with G — c = 1. 



2 DESCRIPTION OF THE GRAVITATIONAL 
WAVE SIGNAL 

The theory of the GW signal produced by the superposition 
of radiation from a large number of individual sources, was 
extensively presented in SVC08 and SVV09, here we review 
the basic concepts, deriving the GW signal for gas driven 
mergers. 

The dimensionless characteristic amplitude of the GW 
background produced by a population of binaries with a 
range of mass es mi and m,2 and redshifts z is given by 
|Phinnevll20"0lh 



h c(f) = / dzd7nidm 2 



d 3 n 



dzdmidm,2 1 + z din f r 



(1) 



where dE gw / 'din f r is the total emitted GW energy per log- 
arithmic frequency interval in the comoving binary rest- 
frame, f r = (1 + z)f is the rest-frame frequency, / is the 
observed frequency, n is the comoving number density of 
sources, and the 1/(1 + z) factor accounts for the redshift 
of the observed GW energy. The characteristic GW ampli- 
tude h c is related to the present day total energy in GWs as 

PGW=fJ7^(/)d/. 



2.1 GW-driven inspirals 

To the leading quadrupole order, for circular purely GW- 
driven binaries orbiting far outside the innermost stable cir- 
cular orbit (ISCO), Eq. (fTJ) can be evaluated assuming that 
i?gw = Spot = —mim,2/a is equal to the Newtonian poten- 
tial energy of the binary, an d that f r is eq ual to twice the 
Keplerian orbital frequency (|Phinney||200ll ) . In this case 



— — - = -/x(7rAf/ r 
d In f r 3 



2/3 1 / r \2/3 . ,5/3 
= jWr) M 



(2) 



for f r < /isco = l/(6 3/2 7rAf). Here /j = mim 2 /M, M = 
mi + m2, and ftA 5 ^ 3 = /xM 2//3 are the reduced, total, and 
chirp masses for a binary, and /isco is the GW frequency 
at ISCO. Substituting into Eq. flTJ, 



hl(f) 



4/ 



-4/3 



3ttV3 



djzdmidm2 



d 3 n 



2/3 



dzdmxdmi (1 + z) 1 / 3 



(3) 



It is also useful to examine the number of M BHBs and their 
respective contributions to the total signal (|Phinnevll200ll ; 
SVC08) . Eq. © can be rewritten as 



K{f) 



dzdmidm2 



d 4 N 



-h\M,Z,fr), (4) 



dm\dm-2dzd In f r 

where iV is the number of sources, which we can calculate 
from a comoving merger rate density as explained below, 
and 

h(M,Z,f r ) = * A *Mfr) 2/3 , 



W 1 / 2 d 



LIZ 



(5) 



is the sky and polarization averaged characteristic GW 
strain amplitude of a single binary with chirp mass M, at 
the particular orbital radius corresponding to f r . 

We generate the distribution d 4 N/(dmidrn2dzdlnfr) 
corresponding to a comoving merger rate densitjij 
d 4 N/ (dmidm2dt r dV c ) (see Sec. I3.1[l . assuming that the 
number of binaries emitting in the interval In f r is propor- 
tional to the time the binary spends at that frequency, 



d 4 N 



d 4 N 



6V C dz dt r 



dm\dm2dzd\n f r dm\dm2dt r dV c dz dt r dln/ r 



(6) 



where dV c /dz and dz/dt r are given by the standard cos- 
mological relations between c omov ing volume, redshift, and 
time, given in e.g. IPhinnevI (j200fh. The last factor can be 
expressed using the residence time t ICS = a(da/dt r )~ 1 the 
binary spends at a particular semimajor axis as 



dt r 




dln/ r 





dt r d In a 



d hi a d In f r 



_ 2 ^ 



(7) 



where Kepler's law, a = M(tvM / r )~ 2//3 , was used to obtain 
dlna/dln/ r = 2/3, and the residence time for a purely GW- 
driven evolution is 

dt r 



tr 



' r 



dlna = U M ^T»». 



(8) 



In summary, the distribution of sources in Eq. Q becomes 
d 4 N 2 d 4 N dV c dz 



dmidm2dzdln f r 3 dm\diri2dt r dV c dz dt 



■tr 



(0) 



1 In practice d 4 N/ (dmidm2dt r dV c ) is a function of mi, rri2, and 
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Equations (0 and (|4]^9]) are equivalent, but ([4H9J) are prac- 
tical to generate discrete Monte Carlo realizations of a 
given source population. Moreover, Equations ((4H2]) provide 
a transparent interpretation of ([3]). The total RMS back- 
ground, h c oc \/~Nh oc /~ 2//3 , comes about because the mean 
number of binaries per frequency bin is N oc tf™ s °c fr~ S ^ 3 
and each binary generates an RMS strain h oc J 2 '' 3 . The 
scaling h c oc /~ 2//3 is a consequence of averaging over the 
local inspiral episodes of merging binaries in the GW driven 
regime, but is completely independent on the overall cos- 
mological merger history or on the involved MBH masses. 
The latter affects only the overall constant of proportional- 
ity. Eq. ([3]) also shows that this scaling constant is insen- 
sitive to the cosmological redshift distribution of mergers, 
h c oc (1 + z) 1//6 , as well as the number of minor mergers 
with n <C M, once the t otal mass in sa tellites merging with 
BHs of mass M is fixed (|Phinnevll200lh . However, the back- 
ground is sensitive to the assembly scenarios of major merg- 
ers (SVV09). More importantly, as shown in SVC08, the 
actual GW signal in any single realization of inspiralling bi- 
naries is qualitatively different from h c oc /~ 2//3 , as a small 
discrete number of individual massive binaries dominate the 
nHz GW-background, creating a very spiky GW spectrum. 
We further discuss the discrete nature of the signal in Sec- 
tion 12.31 below. 



2.2 Gas-driven inspirals 

Let us now derive the GW background for an arbitrary 
model of binary evolution. We can derive the background 
using the residence time t res the binary spends at each semi- 
major axis a, where t les < tf™ in the gas driven phase (cf. 
Eq. [8] for the purely GW driven case). We define t rcs for 
various accretion-disk models in Section [3.21 Generally, the 
emitted GW spectrum is 



dE s 



dE s 



dt r 



2 dE s 



d In f r 



dt r d In f r 



3 dt r 



-trc 



(10) 



The second equality follows the definition of t res and Kepler's 
law (see Eq. [7]). The emitted power d_E gw /dt r depends only 
on the masses and the geometry of the orbit, but is indepen- 
dent of the global migration rate of the binary, and therefore 
it is the same as in the pure GW driven case. The effects of 
migration is fully encoded in t rcs . Plugging back in Eq. |T]), 
the mean square signal in the gas driven phase is 



K(f) 



if 



-4/3 



37T 1 / 3 

t Ies (M,fl,f r ) 
tfZ(MJr) 



dzdmidm,2 



fj,M 



2/3 



dzdm\dm2 (1 + z) 1 / 3 

(11) 



A comparison with ([3} shows that the RMS signal is attenu- 
ated in the gas driven phase by t Ies (M, \x, f r ) /fJS {M , fr)- 
This factor is a complicated function of f r , that behaves dif- 
ferently for different masses and mass ratios of the binaries; 
the overall spectrum is no longer a powerlaw. 

We generate Monte Carlo realizations of the GW signal 
by sampling the population of the inspiralling systems. To 
do this, it is sufficient to recognize that the derivation given 
by Eqs. (j3H§| remains valid in the gas drive phase if using 
the appropriate t las (M , fi, f r ) , since the contribution given 
by each individual source to the signal is the same. The net 



GW spectrum changes because of a reduction in the number 
of sources in the gas driven case. 

Note that the individual GW signal given by Eq. (JSJ) de- 
pends on three parameters only: M, z, and f r - This implies 
that signals from sources with the same observed frequency, 
redshift, and chirp mass, (f r , z,M), but different mass ra- 
tios, q, are totally indistinguishablfl We can make use of 
this property and reduce the number of independent param- 
eters in the distribution by integrating over the mass ratio 
in Eq. © 



dMdzdlnfr 



dq 



dm\dm-2dzd In f r 



<9(mi, 7712) 



8(M,q) 



(12) 



Note, that this step is different for the GW and gas driven 
cases, because d 4 N/(dmidm2dzdln f r ) is proportional to 
trcs in Eq. JjJ. Here \d(mi, m2)/d{M, q)\ is the determi- 
nant of the Jacobian matrix corresponding to the variable 
change from (7711,7712) to (M,q). With Eq. @ and (JT^J we 
derive d 3 N/(dMdzd\nf r ) for any gas driven model given 
by t ICS (M, fi, fr), and draw Monte Carlo samples of inspi- 
ralling binaries from this distribution when generating the 
GW signal. 



2.3 Statistical characterization of the signal 

In observations with PTAs, radio-pulsars are monitored 
weekly for total periods of several years. Assuming a re- 
peated observation in uniform At time intervals for a total 
time T, the maximum and minimum resolvable frequencies 
are /max = 1/(2 At), corresponding to the Nyquist frequency, 
and / m in = 1/T. The observed GW spectrum is therefore 
discretely sampled in bins of A/ = /mm. For circular orbits, 
the frequency of the GWs is twice the orbital frequency. 

Let us examine whether the sources' GW frequency 
evolves during the observation relative to the size of the fre- 
quency bins. Writing the frequency shift during an observa- 
tion time T as A/ evo] « fT = (dln//dlna)(dlna/dt)/T = 
|/T/[(1 + z)t TQS ] and considering the frequency resolution 
bin to be A/bin = 1/T, then the frequency evolution relative 
to the frequency resolution bin is 



A/ c . 



fT 2 



A/ b i n 2 (1 + Z)t r 



^M^f^T?^, (13) 

L ~\~ Ai ^res 



where Als.s is the chirp mass in units of 10 Mq, /50 is 
the frequency in units of 50 nHz, T10 is the observation time 
in unit of 10 years, and for the second equality we have 
used equation ((8]). Equation (|13|) shows that typical binaries 
contribute to a single frequency bin as stationary sources in 
the GW-driven regime0 We shall demonstrate that this is 
also true in the gas driven case (see Figure [1] below). 

We generate iVfc different realizations of the signal (usu- 
ally N k = 1000), i.e. N k realizations of the MBHB popula- 
tion, consistent with Eq. (|9]) (see Section [3]T}. Each of those 
consists of Nb ~ 10 3 — 10 4 binaries producing a relevant 
contribution to the signal, which we label by (Mi, Zi, f ri ), 
i — 1,2, ... , iVb- The total signal (Eq. [4]) in each frequency 



2 This is true only in the angular averaged approximation. We 
neglect the directional sensitivity of PTAs. 



3 The detection of the frequency shift of an M%. 5 = /50 = Z ■ 
source would require an extended observation with TJ>35 yr. 



1 



© 0000 RAS, MNRAS 000, 000-000 



Massive black holes and pulsar timing 5 



resolution bin A/ is evaluated as t he sum of the contri- 
bution s of each individual source fsee lAmaro-Seoane et al.~l 
l)20ld ) for the detailed numerical procedure) 



hl{f) = y^ft-c,» (Mi, Zj,f T 



(14) 



where for each / the sum is over the inspiralling sources 
emitting in the corresponding blue-shifted (i.e. restframe) 
A/ r frequency resolution bin, and h c ,i = hi-s/fiT is the angle 
and polarization averaged GW strain given by Eq. ([5]), mul- 
tiplied by the square root of the number of cycles completed 
in the observation time. For comparison, we also evaluate 
the continuous integral Eq. Q, which represents the RMS 
average of Eq. (|14|l over different realizations of MBHB pop- 
ulations in the limit Nk — > oo. 

Since the mass function of merging binaries is in gen- 
eral quite steep, the relative contributions of the few most 
massive binaries turn out to dominate the background in 
each frequency bin. The total GW signal depends very sen- 
sitively on these rare binaries, and the inferred spectrum is 
very spiky. It is useful to separate the total signal h c into a 
part generated by a population of GW-bright individually 
resolvable sources, and a stochastic level h„, which includes 
the contribution of all the unresolvable, dimmer sources. 
More precisely, in each frequency resolution bin, we find the 
MBHB with the largest hci(Mi, Zi, f ri ), and we define it 
individually resolvable if its signal is stronger than the total 
contribution of all the other sources in that particular fre- 
quency bin. The stochastic level is consequently defined by 
adding up only the unresolvable sources in Eq. (fT4)) . Since 
the signal is dominated by few individual sources in each 
frequency bin, the h 3 (f) distribution obtained over the Nk 
realizations is far from being Gaussian or just even symmet- 
ric. To give an idea of the uncertainty range of h 3 (f), we 
calculate the 10%, 50% and the 90% percentile levels of the 
h s (f) distribution of the Nk different realizations. 

The separation of individually resolvable sources is use- 
ful for several reasons (see SVC; SVV09). First, it is use- 
ful from a statistical point of view for understanding the 
variance of the expected GW spectrum among various re- 
alizations of the inspiralling MBHBs. The discrete nature 
of the resolvable sources allows a different statistical anal- 
ysis than for the smooth background level corresponding 
to the stochastic level, h B (f). The individually resolvable 
signals could also be important observationally. A suffi- 
ciently GW-bright resolvable binary allows to measure the 
GW polarization using PT As, and give information on the 
sky position of the binary (|Sesana which 
might be used to search for direct electromagnetic signa- 
tures like periodically variable AGN activity (HKM09). A 
coincident detection of GWs and electromagnetic emission 
of the same binary would have far reaching consequences 
in fundamental physics, cosmolo gy, and black hole physics 
l|Kocsis. Haiman. fc Menoull20Qg| ). 



sponding to the sky position and polarization averaged delay 
in the time of arrivals of consequent pulses due to GWs, 



8tc(f) = 



ho(f) 
2nf ■ 



(15) 



The pulsar timing residuals expected from an individual 
stationary GW source is derived in section 3 of SVV09 in 
detail. The corresponding measurement can be represented, 
in the time domain, with a residual: 



5t{t) = r(t) + 6t N (t) 



(16) 



where r(t) is the contribution due to the GW source (which 
accumulates continuously with observing time t, see below), 
and St^ (t) represents random fluctuations due to noise. The 
latter is the superposition of the intrinsic noise in the mea- 
surements and the GW stochastic level from the whole pop- 
ulation of MBHBs, with a root-mean-square (RMS) value 



Sti 



s(f) = {stUf)) = K(f) + 5t2 a (f)- 



(17) 



where t p (f) is the RMS instrumental and astrophysical noise 
corresponding to the given pulsar, and St s (f) — h s (f)/(2nf) 
is due to the RMS stochastic GW background of unresolved 
MBHBs, as defined in the previous section. 

The sky angle and orbital orientation-averaged signal- 
to-noise ratio (SNR) at which one MBHB, radiating at 
(GW) frequency /, can be detected using a single pulsar 
with matched filtering is 



SNR 



sjUf) 
stl ims (f) 



(18) 



Here <5i gw (/) is the root-mean-squared timing residual signal 
resulting from GWs emitted by the individual stationary 
source over the observation time T defined as: 



St s-»(f) = 



8 h(M,z,f) 



IT 



(19) 



15 2tt/ 

where h(M, z, /) is the angle and polarization averaged GW 
strain amplitude given by Eq. ([SJ, the prefactor ^/8/15 av- 
erages the observed signal over the "antenna beam pattern" 
of the array (Eq. (21) in SWOQ^), and the yJJT term ac- 
counts for the residual build-up with the number of cycles. 
For N p number of pulsars, the total detection SNR, of an 
individually resolvable MBHB is the RMS of the contribu- 
tions of individual pulsars given by (|18[) . For N p identical 
pulsars, the effective noise level is therefore attenuated by 
N-^\ 

In the following we will represent the overall GW signal 
and the stochastic background by using either their charac- 
teristic amplitudes, h c (f) and h s {f), or the corresponding 
characteristic timing residuals, St c (f) and St s (f), according 
to equation (I15|l . We study the detection significance of indi- 
vidually resolvable sources and the distribution of their num- 
bers as a function of the induced <5t gw . For each Monte Carlo 
realization of the emitting MBHB population, we count the 
cumulative number of all (Nt) and resolvable (N r ) sources 



2.4 Timing Residuals 

In general, the characteristic GW amplitudes (either of a 
stochastic background or of a resolvable source) can be 
translated into the pulsar timing language by converting 
h c (f) into a "characteristic timing residual" St c (f) corre- 



N t /r(8tgw) = 



ON 



t/r 



dSt' 



St' 



Note that that the square root in the prefactor 1 
ing in Eq. (20) of SVV09 because of a typo there. 



(20) 



'8/15 is miss- 
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where the integral is either over all sources or only the indi- 
vidually resolvable sources (i.e. restricted to those that pro- 
duce residuals above the RMS stochastic level, see Sec. I2J3JI . 



3 THE EMITTING BINARY POPULATION 

We calculate the GW signal by generating a catalogue of 
binaries consistent with Eq. (|9]). This requires (i) a model 
for the comoving merger rate density of coalescing MBHBs, 
d 4 N/(dmidm2dt r dV c ), and (ii) a model for the evolution 
of individual inspiralling binaries t les {M, fi, f r ). These two 
items are the subjects of the next two subsections below. 



3.1 Population of coalescing massive black hole 
binaries 

We use the population models described in section 2 of 
SVV09, the reader is deferred to that paper for full de- 
tails. We extract ca talogs of merging bin aries from the semi- 
analytical m odel oflBertone et alj (|2007T ) applied to the Mil- 
lennium run (|Springel et alj|2005l ). We then associate a cen- 
tral MBH to each merging galaxy in our catalogue. We ex- 
plored a total of 12 models, combining four MBH-bulge pre- 
scriptions found in the literature with three different ac- 
cretion scenarios during mergers. The twelve models are 
listed in table 1 of SVV09. In the present study we shall 
use the Tu-SA population model as our default case. In this 
model, the MBH masses in the merging galaxies correlate 
wi th the masses o f the bulges following the relation reported 
in iTundo et al~l (|2007t) . and accretion is efficient onto the 
more massive black hole, before the final coalescence of the 
binary. Since the comparison among different MBHB popu- 
lation models is not the main purpose of this study, we will 
present results only for this model. However, we tested our 
dynamic scenario on other population models presented in 
SVV09, finding no major differences for alternative models. 

Assigning a MBH to each galaxy, we obtain a cata- 
logue of mergers labelled by MBH masses and redshift. From 
this, we generate the merger rate per comoving volume, 
d 4 N/ (dmidm2dt r dV c ) . In practice, due to the large number 
of mergers in the simulation, this is a finely resolved continu- 
ous function, describing the merger rate density as a function 
of mi, m2, and z. After plugging into Eqs. (0 and ([12]) . we 
can obtain the continuous distribution d 3 N/(dMdzdln f r ) 
one would observe in an "ideal snapshot" of the whole sky. 
We then sample this distribution to generate random Monte 
Carlo realizations of the GW signal. In summary, the chosen 
MBHB population model fixes the cosmological merger his- 
tory, i.e. the function d 4 N/(dm 1 dm2dt r dV c ), and "Monte 
Carlo sampling" refers to first choosing a realization of the 
cosmological merger history (i.e. generate Nb number of bi- 
nary masses and redshiftqj) and then assigning an orbital 
frequency (or time to merger) to each binary. In addition to 
our fiducial merger rate, we also examine for the first time 
the situation where minor mergers do not contribute to the 
coalescence rate of the central MBHs. This is motivated by 



recent numerical simulations indicating that minor mergers 
with mass ratios q < 0.1 lead to tidal stripping of the merg- 
ing satellite, and the resulting co re does not sink effici ently 
to the centre of the host galaxy |Callegari et al.ll2009h . We 
identify the mass ratio of merging galaxies using the mass 
of the stellar components, to avoid complications due to the 
tidal stripping of merging dark matter halos. In practice, we 
find that suppressing all the minor mergers does not affect 
the resulting GW signal, implying that the contribution of 
mergers involving dwarf galaxies is negligible. 

3.2 Binary evolution in massive circumbinary 
disks 

We adopt the simple analytical models of HKM09 to de- 
scribe the dynamical evolution of MBHBs in a geomet- 
rically thin circumbinary accretion disk. Here we provide 
only a short summary of the gas-driven models, highlight 
their main assumptions, and defer the reader to section 2 of 
HKM09 and references therein for more details. 

For the typical MBH masses and separations we are con- 
sidering, the tidal torque from the binary dominates over the 
viscous torque in the disk, opening a gap in the gas distri- 
bution. A spiral density wave is excited in the disk which 
torques the binary and pushes it inward. First, when the 
binary separation is relatively large, and the local disk mass 
is larger than the mass of the secondary, then the secondary 
migrates inward with the radial velocity of the accreting gas, 

Cs = _^ = 2 ![ r|Eo 
a v M 

where ro is the radius of the gap (typically ro ~ 1.5a where 
a is the semimajor axis of the binary), Eo is the local sur- 
face density of the disk with no secondary, and M is the 
accretion rate far from the binary (the v index denotes vis- 
cous evolution). This is analogous to disk dominated Type-II 
planetary migration. More typically, however, the binary is 
more massive than the local disk mass. In this case, analo- 
gous to secondary-dominated Type-II migration, the angular 
momentum of the binary is absorbed less efficiently by the 
gas outside the gap, and the evolution slows down according 
to i|Sver fc Clarke! 1 1995T ) 

d° = -A=fe fel t„ (22) 
asc 

where gs is a measure of the lack of local disk mass domi- 
nance 



IB 



i 'is • 



(23) 



5 Here TVj, is chosen randomly for each distribution, it can vary for 
different realizations of the population according to a Gaussian 
distribution with a = 1/ \f (N^) around the mean (iVj,). 



which is less than unity in this case; fi is the reduced mass 
of the binary, and k\ is a constant that depends on the 
surface density - accretion rate powerlaw index. Here ki = 
if qs > 1; otherwise, k\ = 7/17 at large separations, if 
the disk opacity at the gap boundary is dominated by free- 
free absorption, and ki = 3/8 closer in, if the opacity is 
dominated by electron scattering. 

In the secondary-dominated regime, eq. (|22p assumes 
that the migration is very slow, so that the binary can 
influence the surface density very far upstream, and an 
approximate steady state is reached where the gas den- 
sity is enhanced outside the gap by a fa ctor q„ k rela- 
tive to the gas density with no secondary (|Sver fc Clarkel 
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1 19951 ). To examine the sensitivity of our conclusions to these 
as sumptions, we consider an addi t ional nonsteady model 
of llvanov. Papaloizou. fc PolnarevI (1 19991 ). In their model, 
which is also axisymmetric by construction, they assumed 
that the binary torques are concentrated in a narrow ring 
near the edge of the cavity, which results in a time-dependent 
pile up of material near the cavity. The migration slows down 



k 2 



2M 



\aoJ 



1/2 



1 + 5(1- ^fafc 



(24) 



where M is the constant accretion rate far outside the bi- 
nary, ao is the initial semimajor axis where the local disk 
mass is just equal to the secondary (i.e. qs — 1), 5 and fe 
are constants which depend on the density-accretion rate 
profile of the disk (S = 4.1 and 6.1, k 2 = 0.29 and 0.26 for 
free- free and electron scattering opacity, respectively). Note 
that as the separation decreases well below its initial value, 
a <C fflo, the curly bracket is a-independent, implying that 



^res OC Q ^ 



OC t 



Since the migration rate is proportional to the gas sur- 
face density and accretion rate, it is sensitive to th e struc - 
ture of the accretion disk. Following ISver fc Clarke! l|l995h . 
we estimate Eo with that of a steady accretion disk of a 
single accreting BH, but consider two different models for 
that, a an d /3-disks (see HKM0 9 for e xplicit formulae). For 
the classic IShakura fc Sunvaevl (|l973l ) a-disk, the viscosity 
is proportional to the total (gas+radiation) pressure of the 
disk. Until very recently, this model, if radiation pressure 
dominate d, has been thought to be thermally an d viscously 
unstable (|Lightman fc Eardlevlll974l ; |Piranlll978l ). In the al- 
ternative /3-model, the viscosity is proportional to the gas 
pressure onljQ, and it is stable in both sense. The nature of 
viscosity is not well understood to predict which of these pre- 
scriptions lies closer to reali ty. Recent numerical magneto - 
hydrodynamic simulations l|Hirose. Krolik. fc Blaes 1 120091 ) 
suggests that the thermal instability is avoided in radia- 
tion pressure dominated situations because stress fluctua- 
tions lead the associated pressure fluctuations, and seem 
to favor the a prescription over the /3-model. We carry 
out all calculations for both models, but consider the a 
prescription as our fiducial disk model. In both cases, the 
model is uniquely determined by three parameters: the cen- 
tral BH mass, the accretion rate far from the binary M, 
and the a viscosity parameter. The exact value of these pa- 
rameters is not well known. Observations of luminous AGN 
imply an accretion rate around rh = M / Msdd = (0.1-1) 
with a statistical increase towards higher qu asar luminosities 
(|Kollmeier et all 120061 ; iTrump et alj |2009| ). Here M Edd = 
iEdd/(??c 2 ) is the Eddington accretion rate for rj = 10% ra- 
diative efficiency, where I/Edd is the Eddington luminosity. 
Observations of outbursts in binaries with an accreting white 
dwarf, neutron star, or stella r black hole imply a = 0.2-0.4 
jDubus. Hameurv. fc Lasotal l200ll ; iKing. Pringle. fc Livid 
12001 and references therein). Theoretical limits based on 
simulations of magneto-hydrodynamic turbulence around 



black holes are i nconclusive, but are consistent with a in 
the range 0.01-1 (|Pessah. Chan, fc Psaltisll2007l ). It is how- 
ever unclear whether these numbers are directly applica- 
ble to circumbinary MBH systems. The binary exerts a 
torque that pushes the gas away on averag e, and conse- 
quen t ly might reduce the accre tion rate (see ILubow et al.l 
Il999t ILubow fc D'Angeld l200d. in the planetary context , 



and IMacFadven fc Milosavlievicl 120081 : ICuadra et ail 120091 . 
for MBH binaries). We explore several choices covering all 
6 combinations with rh = {0.1, 0.3} and a = {0.01, 0.1, 0.3} 
for both a and /3-disks, respectively. Motivated by the con- 
siderations outlined above, we highlight the a-disk with 
rh = 0.3 and a = 0.3 as our default model. 

We regard tf2. an d tJ-c S P as lower and upper lim- 
its of the true residence time during secondary domi- 
nated gas driven migration. However, we caution that 
these estimates are subject to many uncertainties re- 
lated to the complexity of binary accretion disk and mi- 
gration physi cs. Other mode l s for the binary-dis k in- 
teraction (i.e. iHavasakil 120091 ; iHavasaki et all |201Ch lead 
to a slower migration rate. Two dimensional hydrody- 
namic simulations show that the flow into the gap is 



nary orbit becomes eccentric ( 


Lubow & D'Angeld 120061: 


IMacFadven & Milosavlievicl 2008 


: ICuadra et all 20091). The 



accretion rate and inspiral velocity in reality may be higher 
than in the axisymmetric approximation. Resonant interac- 
tions may lead an enhanced angular momentum transport 
jGoldreich fc Tremainlll979l ). For unequal masses or thick 
disks, there may be significant inflow into the gap. In this 
case, corotation torques need to be considered, and the mi- 
gration is expected to be much faster inwar d or outward, 
analo gous to Type I planetary migration (|Tanaka et al.l 
|2002| ). Disk thermodyn amics may also significantly influ- 
ence the migration rate (Paardekooper fc Papaloizou |[2009l ; 
iPaardekooper et alj|2010h . 

These steady state accretion disk models also assume 
that the self-gravity of the disk is negligible, the disk is op- 
tically thick, and the temperature of the gas is large enough 
(TJ>10 4 K) that the opacity corresponds to that of ionized 
gas. These assumptions break down at radial distances cor- 
responding to orbital periods of around t or b^(5-10) yr for 
BHs with masses M^(10 8 -10 9 ) M Q (see Fig. 1 of HKM09). 
Thus, these disk models are self-consistent for the relatively 
short period binaries, but become suspect for binaries emit- 
ting at the low frequency edge of the PTA window. 

Another shortcoming of these gas-driven migration 
models, is that they were derived assuming the disk sur- 
face density is a decreasing function of radius. This affects 
the specific values of fci and k 2 in eqs. (|22p and (|24|l . This 
assumption is safe for /3-disks in general, as well as for 
a-disks if the gas pressure dominates over radiation pres- 
sure, but it is violated for radiation pressure dominated a— 
disks. However, gas-driven models with binary separation in 
the PTA frequency band and large masses M > 10 8 M 
are in fact radiation pressure dominated near the gap 
It is unclear how k\ and k 2 change for radiation pres- 
sure dominated a-disks. As an approximation, following 



6 The name comes from the definition of viscosity v oc opg aB pJ ot ' 3 



where /3 = 1 for the /3— model, while f) 
both cases a is a free model parameter. 



for the a model. In 



7 Note that the most massive binaries, however, are GW driven, 
and arc not sensitive to the accretion disk. 
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Figure 1. Residence time, t TCB , as a function of the binary pe- 
riod. Different panels correspond to selected model parameters 
as labelled. In each panel the thick curves represent the resi- 
dence time for the gas driven dynamics according to HKM09. 
Solid, long— dashed and short-dashed curves are for mass ratios 
q = 1,0.1,0.01 respectively. The two thin dotted vertical lines 
approximately bracket the PTA observable window. The extrap- 
olated pure GW-driven evolution (tf™ s oc t S is shown as thin 
lines, for comparison. The ratio t r cs/if<S gives the relative de- 
crease in the number of binaries in the gas dominated case com- 
pared to the GW driven case. 



lArmitage Sz Nataraian I (|2005h and HKM09, we extrapolate 
the values of the gas pressure dominated outer regions. 

Substituting the surface density profiles So(r) of a and 
/3-disks in eqs. (|21II22|I . we obtain the inspiral rate of bi- 
naries as they evolve from large to small radii in a gaseous 
medium. Figure [1] shows the corresponding residence times 
for selected MBH masses and disk models. The binary mi- 
gration history, can be divided in three phases. 



(i) Initially, migration occurs on the viscous timescale 



t' ' , 



,5/6 



(ii) Then as the disk mass decreases below the secondary 
mass, the gas-driven migration slows down according to tf^ 
or ires' '• The migration rate in this case is nonuniform, it 
changes as the local accretion disk structure varies due to 
changes in the dominant source of opacity (free-free for 
large separation to electron-scattering for smaller separa- 
tions) and the source of pressure (thermal gas pressure at 
large a to radiation pressure at small a). In particular for 

25/51 ,7/12 



orb 



orb 



the steady ISver fc Clarkel (| 19951 ) models, « 1 
with free-free to electron scattering opacity for thermal pres 
sure support, and oc t 



35/24 .7/12 , 



for radiation pres- 



sure support for a to /3 d isks, respectively; while for the 
time-d ependent models of llvanov. Papaloizou. k, Polnarevl 
the accretion rate approaches t]^ oc t^ 3 . Note that 
the residence time in this regime is generally smaller for a- 
disks than for /3-disks. This is explained by the fact that 
the overall surface density is smaller for radiation pressure 



dominated a-disks, and so to achieve the same net accretion 
rate, the radial gas inflow velocity is larger for a-disks. The 
migration rate is a monotonic function of the radial inflow 
velocity (see eqs. l21l and l22[) . so that the binary is pushed in 
more quickly for a-disks. 

(iii) Finally, when the separation is sufficiently reduced, 
the emission of gravitational waves becomes very efficient 
and determines the inspiral rate according to tf™ B = a/d gw oc 

b ( see e iM- 



4 ,8/3 

a oc t 



Note that for circular orbits, the GW frequency (in 
the binary restframe) is simply f r = 2/r or b. The number 
of binaries at any given t or b is proportional to the resi- 
dence time t res . Therefore, the decrease in t les in the gas 
driven regime compared to the GW driven case implies a 
decrease in the population of MBHs, which ultimately leads 
to the attenuation of the low frequency end of the observ- 
able total GW spectrum. The RMS GW spectrum averaged 
over the whole population of binary inspiral episodes is no 
longer a powerlaw. It is interesting to examine the contribu- 
tions of various evolutionary phases according to Eq. (| 1 1 p . 
h c oc / _2//3 \J t ICS /tf™ s . If all binaries were in the GW driven 
phase then hf oc /" 2/3 . If all were in the secondary- 
dominated Type-II migration regime with a steady radia- 



tion pressure dominated disk 



oc /" 1/16 -/ 3/8 for a and 



/3 disks, respectively, further out oc J 3 / 8 -/ 43 / 102 j n the 
gas pressure dominated regime (with electron scattering ver- 
sus free-free opacity), while h lpp oc / 1//2 asymptotically for 
the nonsteady models, and finally h v c oc / 1//4 in the disk- 
dominated Type-II migration regime. In any case, the GW 
spectrum h c (f) is much shallower in the gas driven phase. 
Note the gas driven phase contributes a nearly flat or an 
increasing spectrum h c (f), very different from the nominal 
/- 2 / 3 GW driven case. In general, the orbital separation 
for more massive objects at a given f r is smaller in terms of 
Schwarzschild radii. Since the transition to gas driven mi- 
gration, when expressed in Schwarzschild radii, is roughly 
independent of mass, it follows then at any given frequency 
bin the more massive objects are typically GW driven and 
lighter objects are gas driven. The total average spectmuQ 
assuming only wet mergerf"!. is between /~ 2/3 and f 1/2 de- 
pending on the ratio of GW driven to gas driven binaries. 
Note however, that the individually resolvable sources are 
typically very massive and are in the GW-driven phase for 
the relevant range of binary separations, and therefore their 
properties should not be modified by gas effects. 



8 Most of the gas driven binaries contributing to the background 
at large frequencies are in this regime. 

9 if averaging over each binary episode but not over the cosmo- 
logical merger tree 

10 Here "wet" refers to gas rich mergers where the dynamics of 
the binary is driven by both the circumbinary disk and GW 
emission depending on the binary separation, while "dry" refers 
to the case where the binary dynamics is driven by GW emission 
only at all radii. Note that in both cases we neglect interactions 
with stars. 
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Figure 2. Components of the GW signal from a population of inspiralling MBHBs. In the left panel we consider all binaries embedded 
in a gaseous «-disk for our default model (all mergers wet), while in the right panel all binaries are purely GW driven (all mergers dry). 
In each panel, the smooth solid line is the RMS total characteristic GW strain h c using the integral expression BH9t (which correspond 
to an average over iVj, — > oo Monte Carlo realizations), the two dashed black lines represent the RMS total signal (upper) and the 
RMS stochastic background level (lower) averaged over TVj, = 1000 Monte Carlo realizations respectively. The jagged blue line displays 
a random Monte Carlo realization of the GW signal. The small black and red triangles show the contributions of the brightest and 
resolvable sources in each frequency bin respectively. The jagged red line is the stochastic GW background for this realization, i.e., once 
the resolvable sources in each frequency bin are subtracted. The green dots label all the systems producing an RMS residual t gw > 0.3 ns 
over T = 10 years. The dotted diagonal lines shows constant t gw levels as a function of frequency. An observation time of 10 years is 
assumed. 



4 RESULTS 

4.1 Description of the signal 

As stated in Sec. 2.3, the relevant frequency band for pulsar 
timing observations, assuming a temporal observation base- 
line T and a time interval between subsequent observations 
At, is between / m i n « 1/T and / max ~ 1/(2 At) with a res- 
olution A/ fts 1/T. In our calculations we assume a default 
duration of T = 10 yr for the PTA campaign with At ~ 1 
week. This gives / m i n Ri 3 x 1CT 9 Hz, / max ~ 10~ 6 Hz and 
A/ « 1/T « 3 x l(T 9 Hz. The simulated signal is com- 
puted by doing a Monte Carlo sampling of the distribution 
d 3 N / (dzdMdln f r ), and adding the GW contribution of 
each individual source. In each frequency bin, we identify 
the individually resolvable sources and the stochastic back- 
ground. We repeat this exercise Nk = 1000 times for the 12 
steady disk models denned in Sec. 13.21 and for the purely 
GW-driven case. In addition, we also calculate the GW sig- 
nal using the integral expression (fJHSl) using the continuous 
distribution function, which corresponds to the RMS aver- 
age of the GW signal in the Nk oo limit. 

The GW signal and its most important ingredients are 
plotted in Figure [2] The left panel shows results for gas- 
driven inspirals using our default a-disk (i.e. "gas on" - all 
mergers wet), the right panel shows results for purely GW- 
driven inspirals for comparison (i.e. "gas off" - all merg- 
ers dry) for the same underlying cosmological MBHB coa- 



lescence rate. A randomly selected Monte Carlo realization 
of the signal is depicted as a dotted blue jagged line. The 
green dots represent the contributions of individual bina- 
ries; systems producing 5t gw (f) > 0.3 ns timing residual are 
shown. In each frequency bin, the brightest source is marked 
by a black triangle. The individually resolvable sources are 
marked by superposed red triangles, and the stochastic level 
from all unresolvable sources is shown with a solid red jagged 
line. Clearly, the GW signal for any single realization is far 
from being smooth. The noisy nature of the signal is due to 
rare massive binaries rising well above the stochastic level. 
Solid black curves in Figure [2] show h c (f), the RMS value 
of the GW signal averaged over the merger episodes, using 
the integral expression (fJfO- The upper black dashed curve 
is the GW signal, averaged over Nk — 1000 Monte Carlo re- 
alizations of the merging systems. This is still noisy, due to 
the finite number of realizations used, but is consistent with 
the integrated average shown by the solid black curve. The 
lower black dashed curve marks the RMS stochastic compo- 
nent, /i s ,rms(/), averaged over the same Nk = 1000 realiza- 
tions. Figure [2] shows that the gaseous disk greatly reduces 
the number of binaries at small frequencies compared to the 
GW-only case, as gas drives the binaries in quickly towards 
the final coalescence. Consequently, the signal is more spiky 
in the presence of gas. There is a clear flattening of the 
RMS spectrum at low frequencies (/ < 1/yr) in the gas- 
driven case compared to the purely GW-driven case. The 
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stochastic level is more suppressed in the wet-only case and 
is much steeper than f~ 2 ^ 3 at high frequencies. Despite the 
spiky signal in a given Monte Carlo realization, the over- 
all shape of the background is well recognizable in the GW 
driven model, while the characterization of the global shape 
of the signal in the gas driven case appears to be less viable. 

Gas driven migration becomes more and more promi- 
nent at large binary separations, corresponding to large or- 
bital times or small GW frequencies. Therefore, to check 
the ultimate maximum impact of gas effects on future PTA 
detections, we simulated the spectrum for a hypothetical 
T = 100 yr observation baseline. Figure [3] shows a realiza- 
tion of the spectrum for such an extended observation, as- 
suming our default gas model (i.e. all mergers wet). Given 
the extended temporal baseline, the minimum observable 
frequency is pushed down to ~ 3 x 10 -10 Hz, where the spec- 
trum of gas driven mergers is considerably flatter. Moreover, 
the frequency resolution bin is then narrower, the number 
of sources per bin is much smaller, and more sources be- 
come resolvable. At the smallest frequencies, the induced 
timing residual can be higher than l/j,s and more than 50 
sources will be individually resolvable at Ins precision level, 
some of them with an SNR as high as 100. These numbers 
are not severely modified if considering purely GW-driven 
dry mergers only (see Fig. [5] below). Interestingly, due to 
the high frequency resolution of such an extended observa- 
tion, the GW frequency of some of the resolvable binaries 
may evolve significantly during the observation (e.g. a bi- 
nary with masses mi = ni2 — 5 x 10 8 Mq or higher are 
nonstationary relative to the bin size at frequencies above 
/>40nHz, see equation [13]). Therefore it may be possi- 
ble to detect the frequency evolution of individually resolv- 
able binaries during such an extended monitoring campaign. 
However, this computation is idealized: it is questionable if 
millisecond pulsars can maintain a ns timing stability over 
such a long timescale, nonetheless, it points out the enor- 
mous capabilities of long term PTA campaigns. 



4.2 Individually resolvable sources 

Let us now examine the prospects for detecting individual 
sources using PTAs. How many sources are expected to be 
individually resolvable? How significant is their detection? 

As explained in Sec. 12. 41 the signal from a specific binary 
can be detected using a PTA if the corresponding timing 
residual <5ig W is above the RMS noise level St^ rms character- 
izing the PTA given by Eq. (|18|) . We should notice, however, 
that this estimate of the number of individually resolvable 
sources is conservative and is likely to provide only a lower 
limit for the following reasons. Firstly, we average over the 
sky position of the binaries and pulsars, while in reality, it 
may be possible to take advantage of the different GW polar- 
ization and amplitude generated by sources in different sky 
positions to deconvolve their signal (even if they have similar 
strength and frequency). Secondly, the brightest source iden- 
tification algorithm can be implemented recursively, after an 
accurate subtraction of the identified sources. However, we 
find that, especially at low frequencies, the distribution of 
GW source amplitudes for various binaries in a single fre- 
quency bin is not strongly hierarchical, so that a recursive 
brightest source finding algorithm shouldn't increase signif- 
icantly the number of resolvable systems. We present here 
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Figure 3. Same as left panel of Figure [2] but for an observation 
lasting 100 years (Averages over the 1000 realizations are not 
shown in this case for clarity). 



results both in terms of the total number (N t ) and resolvable 
systems (N T , see Eg. I2TJ1). 

In Figure [4] we plot the distribution of the number of 
sources (total and resolvable) as a function of timing resid- 
ual, detection frequency, redshift, and chirp mass, found in 
two particular realizations of our default a-disk (all merg- 
ers wet) and in the purely GW-driven models (all mergers 
dry). The figure shows that even though there are much 
more sources in the purely GW-driven case, the number of 
sources rising above the stochastic level is almost the same in 
the purely dry and wet cases. Figure 0] also shows the chirp 
mass, redshift, and frequency distribution for sources above 
1 ns timing level. As was previously shown in SVV09, the 
bulk of the sources are cosmologically nearby (z<,l) with 
masses peaking around M ~ 2 x 10 s Mq . Figure U shows 
that gas dynamics does not introduce a major systematic 
change in the shape of the redshift and the chirp mass dis- 
tributions. The lower left panel shows that gas removes a 
systematically larger fraction of sources at small frequen- 
cies. 

Figure [5] shows the cumulative number of binaries (to- 
tal and resolvable) as a function of timing residual. The up- 
per panel shows results for the a-disk models with m = 0.3 
and different a. The statistics of resolvable sources is almost 
unaffected by the large suppression of the total number of 
sources at a fixed timing residual. For example, in all of our 
models, we expect ~ 2 resolvable sources at a timing level 
of <5t gw = 10 ns, even though the total number of sources 
contributing to the signal at that level spans about an order 
of magnitude among the different models (~ 5 for a — 0.3 
to ~ 50 for binaries driven by GW only). The same is true 
for /3-disk models (central panel) and for nonsteady models 
(lower panel), even though in these cases the total number 
of sources at a particular <5t gw is not reduced dramatically 
by gas effects. This result can be understood with a closer 
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Figure 4. Top-left panel: characteristic amplitude of the timing residuals 54 gw (equation d 1911 ") as a function of frequency; the dots 
are the residuals generated by individual sources and the solid line is the estimated stochastic level of the GW signal. Top-right panel: 
distribution of the number of total (dotted lines) and resolvable (solid lines) sources per logarithmic redshift interval as a function of 
redshift, generating a <5t gw > Ins. Bottom-left panel: distribution of the number of total (dotted lines) and resolvable (solid lines) sources 
per logarithmic frequency interval as a function of the GW frequency, generating a <5t gw > Ins. Bottom-right panel: distribution of the 
total (dotted lines) and individually resolvable (solid lines) number of sources per logarithmic chirp mass interval as a function of chirp 
mass, generating a 8t gvJ > Ins. All the black elements refer to our default disk model, the red elements are for a GW driven MBHB 
population. Distributions are averaged over = 1000 realizations of the MBHB population. 



inspection of Figure [T] Let us focus on the a-disk model. As 
explained in Sec. 3.2, the impact of gas driven dynamics in 
the PTA window is more significant for lower binary masses 
and unequal mass ratios. These are the binaries that build- 
up the bulk of the signal, and its stochastic level is conse- 
quently greatly reduced in the gas driven case. On the other 
hand, the population of high equal mass binaries, which con- 
stitute most of the individually resolvable sources, is almost 
unaffected by the presence of the circumbinary disk, as they 
are already in the GW-driven regime in the relevant range 
of binary periods. 

Figure [5] shows expectations on the detection signifi- 
cance of resolvable sources. The SNR distribution of resolv- 
able sources (equation [18]) are shown as thin lines, account- 
ing for the astrophysical GW noise from the unresolved bi- 
naries, but neglecting the intrinsic noise of the array (i.e. as- 
suming an ideal detector with infinite sensitivity, 5tp(/) = 0, 



in equation I17p . This calculation represents an upper limit 
of the SNR. Thick lines, in Figure [6] plot the SNR consider- 
ing a total detector noise of Ins (appropriate for SKA). The 
figure shows that the expected detection significance of re- 
solvable sources is systematical higher for gas driven models 
with larger a. In general, in an array with Ins sensitivity, we 
may expect a couple of individually resolvable sources with 
SNR > 5. For near future instruments with much worse 
sensitivity, the identification of resolvable sources might be 
more challenging. 

4.3 Stochastic background 

Figure shows the RMS stochastic level (i.e. after sub- 
tracting off the individually resolvable sources, see Sec. 12. 3[) 
for selected steady state gas disk models, averaged over 
Nk = 1000 realizations. The top left panel highlights the 
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Figure 5. Cumulative number of total (Nt(5t gw ), thin lines) and 
individually resolvable (N r (St gvi ), thick lines) sources emitting 
above a given 5f gw threshold as a function of <5f gw . Upper panel: 
a-disk model with rh = 0.3. Lines are for a =0.3 (solid), 0.1 
(long— dashed) and 0.01 (short-dashed). Central panel: same as 
upper panel but for a /3-disk model. Lower panel: nonsteady IPP 
model, assuming a = 0.3 and rh = 0.3 (solid) and 0.1 (long- 
dashed). In all the panels, the dotted lines refer, for comparison, to 
the GW driven model. All the distributions refer to the ensemble 
mean computed over all JVj, = 1000 realizations of the MBHB 
population. 
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Figure 6. SNR distribution of individually resolvable sources. 
Thin and thick curves refer to neglecting the instrumental noise 
or using a St gvr = 1 ns timing precision, respectively. Linestyle as 
in the upper panel of Figure \E\ 
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Figure 7. Influence of the gas driven dynamics on h c (thin lines) 
and on h s (thick lines). Upper left panel: GW driven dynamics 
versus gas driven dynamics for our default disk model. Upper right 
panel: a = 0.3 (solid), 0.1 (long-dashed), 0.01 (short-dashed), for 
an a-disk with rh = 0.3. Lower left panel: same as the upper right 
panel, considering an a— disk with rh = 0.1. Lower right panel: 
same as the upper right panel, considering a /3-disk with m = 0.3. 
The two dotted lines in each panel represent the sensitivity of the 
complete PPTA survey and an indicative sensitivity of Ins for the 
SKA. 



difference between wet and dry models, the top right and 
bottom left panels collect different a-disk models and the 
lower right panel is for selected /3-disk models (nonsteady 
models, not shown here, give results basically identical to 
/3-disk models). The different line styles show the effect of 
changing the a parameter of the disk. We also show the RMS 
total signal level, which is exactly proportional to f~ 2 ^ 3 in 
the dry case. In general, the stochastic level matches the 
total signal level at low frequencies, but is increasingly sup- 
pressed for frequencies above /i>10 -8 Hz for all of our mod- 
els. At sufficiently large frequencies GW emission dominates 
even for wet mergers, and both the RMS total signal and 
the stochastic level approaches the purely GW-driven case. 
However, at small frequencies, a significant fraction of bina- 
ries is driven by gas, and the signal is attenuated and the 
spectrum is less steep compared to the dry case. Figure [7] 
shows that gas-driven migration suppresses the stochastic 
background significantly, by a factor of 5 for our standard 
disk model below 10 -8 Hz. The suppression of the total and 
stochastic levels is a strong function of the model param- 
eters. Interestingly, there is almost no suppression for j3- 
disks, or for a-disks with a small accretion rate and/or a 
small a value. In these cases, the local disk mass is smaller, 
resulting in longer viscous timescales, and hence the popula- 
tion of widely separated binaries is not reduced significantly. 

Figure [8] quantifies how the stochastic background 
changes among different Monte Carlo realizations, showing 
the range attained by 10-90% of all Nt = 1000 realizations. 
The variance is not the product of uncertainties related to 
the parameters of the adopted cosmological and dynami- 
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Figure 8. Variance of the expected stochastic level of the signal 
as a function of a. In all the panels we assume the a— disk model 
with rh = 0.3 (except for the lower right panel, referring to the 
GW driven model). Solid lines represent the median h B over 1000 
Monte Carlo realizations, while the shaded area enclosed within 
the two dashed lines is the 10%— 90% confidence region . The two 
dotted lines in each panel represent the sensitivity of the complete 
PPTA survey and an indicative sensitivity of Ins for the SKA. 

cal model, but is purely determined by the small number 
statistics of sources per frequency bin, intrinsic to the source 
distribution. At f — 10 _8 Hz, the variance of the signal pro- 
duced by our default disk model is ~ 0.5dex. Decreasing a to 
0.1 and 0.01, the variance drops to ~ 0.35dex and ~ 0.25dex 
respectively. In the case of all mergers driven by GWs, the 
variance at the same frequency is only ~ 0.15dex. This 
means that, contrary to the GW driven model (SVC08), 
it is impossible to predict the stochastic level of the signal 
accurately for our default gas driven model. For a linear 
fit, any power law in the range /~ 0,3 — / — 1,1 is acceptable 
within the level of variance of the signal in the frequency 
range 3 — 30 nHz . 



5 DISCUSSION AND CONCLUSIONS 

The presence of a strong nHz gravitational wave signal from 
a cosmological population of massive black hole binaries, is 
a clear prediction of hierarchical models of structure forma- 
tion, where galaxy evolution proceeds through a sequence 
of merger events. The detailed nature of the signal depends 
however on a number of uncertain factors: the MBHB mass 
function, the cosmological merger rate, the detailed evolu- 
tion of binaries, and so on. In particular, MBHB dynamics 
determines the number of sources emitting at any given fre- 
quency, and thus the overall shape and strength of the signal. 
Previous works on the subject (e.g. SVC08 and SVV09) con- 
sidered the case of GW driven binaries only. In this paper we 
studied the impact of gas driven massive black hole binary 
dynamics on the nHz gravitational wave signal detectable 
with pulsar timing arrays. This is relevant because in any 



merger event, cold gas is efficiently funnelled toward the cen- 
tre of the merger remnant, providing a large supply of gas to 
the MBHB formed following the galaxy interaction. Even a 
percent of the galaxy mass in cold gas funnelled toward the 
centre, is much larger than the masses of the putative MBHs 
involved in the merger so that the newly-formed binary evo- 
lution may be driven by gas until few thousand years before 
final coalescence. 

To conduct our study, we coupled models for gas driven 
inspirals of HKM09 to the MBHB population models pre- 
sented in SVV09. Our simulations cover a large variety of 
steady-state and quasistationary one-zone disk models for 
the MBHB-disk dynamical interactions (with an extensive 
exploration of the a viscosity parameter and rh for a and 
/3-disks), along with several different prescriptions for the 
merging MBHB population (four different black hole mass- 
galaxy bulge relations coupled with three different accre- 
tion recipes). The differences with respect to the purely GW 
driven models (presented in SVV09) are qualitatively similar 
for all the considered MBHB populations, we thus presented 
the results for the Tu-SA population model only, focusing on 
the impact of the different disk models. Our main findings 
can be summarized as follows: 

• The effect of gas driven dynamics may or may not 
be important depending on the properties of the cir- 
cumbinary disk. A robust result is that if the viscosity 
is proportional to the gas pressure only (/3-disk models), 
there is basically no effect on the GW signal, indepen- 
dently of the other disk parameters. Similarly, there is 
a tiny effect for time-dependent m igration models (i.e. 
Ilvanov. Papaloizou. fc Polnarevl Il999l ) , where gas piles up 
and the accretion rate decreases as the binary hard- 
ens. However, if the viscosity is pro portional to the to- 
tal ( gas+radiation) pressure (a-disk, IShakura fc Sunvaevl 
Il973h . then the GW signal can be significantly af- 
fected for certain st eady state circumbinary disk models 
(|Sver fc Clarke! 1 19951 ) with a>0.1 and m>0.3. This differ- 
ence is explained by the fact that the gas-driven inspiral 
rate is determined by the radial gas inflow velocity, which is 
significantly faster for radiation pressure dominated a-disks 
for a fixed accretion rate, as they are more dilute. There- 
fore gas effects can dominate over the GW-driven inspiral 
rate for a disks at relatively small binary separations corre- 
sponding to the PTA frequency band. 

• With respect to the GW driven case, the presence of 
massive circumbinary disks affects the population of low- 
unequal mass binaries predominantly (M < 10 8 M© , q < 
0.1), causing a significant suppression of the stochastic level 
of the signal, but leaving the number and strength of mas- 
sive individually resolvable sources basically unaffected. In 
our default model (a = 0.3, rh — 0.3), the stochastic back- 
ground is suppressed by a factor of ~ 5 at / < 10 _8 Hz. 
This suppression factor decreases for smaller rh and a. The 
stochastic level is not suppressed significantly for nonsteady 
disks, for /3-disks (arbitrary rh and a) , and for a-disks with 
either m<J0.1 and arbitrary a, or a^O.Ol with arbitrary 
rh. About 10 individual sources are resolvable at 1 ns tim- 
ing level, independently of the adopted disk model. 

• All the results shown here for the Tu-SA model, hold 
for every other MBHB population model we tested. There 
is a certain level of degeneracy between disk dynamics and 
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MBHB mass function: the stochastic level given by a popu- 
lation of heavy binaries evolving by gas dynamics, can mimic 
that of a population of lighter binaries that are driven by 
GWs only. However, the variance of the signal would be 
much bigger in the former case, because the signal is pro- 
duced by fewer massive sources. 

• The detection of GWs emitted by MBHBs embedded in 
gaseous disks with high viscosity and accretion rate, may be 
very challenging for relatively short term PTA campaigns 
like the PPTA. In fact, we find that most of the 12 MBHB 
population models tested in SVV09 would not produce a 
stochastic signal detectable by the PPTA (i.e. the signal is 
a factor of three below the PPTA capabilities for the Tu-SA 
model). However, long term projects like the SKA, which 
aim to nanosecond sensitivities, are expected to be able 
to detect the GW signal, resolving a handful of individual 
sources with high significance. 

A word of caution should be spent to stress the fact 
that our models are idealized in many ways. We considered 
circular binaries only. If most systems were significantly ec- 
centric, then the overall signal would be modified by the 
multi-harmonic emission of each individual source. More- 
over, we only considered radiatively efficient, geometrically 
thin, one-zone, steady-state and quasistationary accretion 
disks. The most massive but gas driven binaries, emitting at 
low GW frequencies / gw ^5nHz, are radiation pressure dom- 
inated, but the migration rate estimates had to be extrapo- 
lated using the scaling exponents for gas pressure dominated 
disks. For even wider massive binaries (/ gw near the low fre- 
quency observation limit) the disks are marginally Toomre 
unstable. We used simple models for the binary disk inter- 
action, scaling the Type-II planetary migration formulae to 
MBHBs. The re sults are sensitive to t he models: the steady- 
state models of ISver fc Clarke] l|l995h lead to a decrease in 
the stochastic nHz backgro und, while the more sophisticated 
time-d ependent mod els of Ivanov. Papaloizou. fc Polnarevl 
(|l999l ) as well as the jHavasaki 1120091 : lHavasaki et al.ll2oTof ) 
models lead to almost no effect (see also Sec. 13.21 for a list 
of caveats). Further studies should examine the accuracy of 
these approximations for gas driven migration in circumbi- 
nary disks around MBHBs. Finally, it is likely that not all 
the binaries are gas driven on their way to the coalescence, 
and the efficiency of the disk-binary coupling may vary from 
merger to merger depending on the environmental condi- 
tions, which may modify the properties of the expected sig- 
nal as well. Nonetheless, our calculations provide clear pre- 
dictions for the possible attenuation of the stochastic GW 
background, which may be confirmed or discarded by ongo- 
ing and forthcoming pulsar timing arrays. 
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